# step7 Scaling the data
# all.genes <- rownames(pbmc)
# pbmc <- ScaleData(pbmc, features = all.genes)


(1) rownames()
这是一个S4风格的函数，解析10讲过了，就是 实现一个 dimnames()，就相当于实现了 rownames() 和 colnames() 2个函数。


./seurat-object-4.0.4/R/seurat.R:2197:#' rownames(pbmc_small)

dimnames.Seurat <- function(x) {
  x <- UpdateSlots(object = x)
  return(dimnames(x = x[[DefaultAssay(object = x)]]))
}




(2) ScaleData()
$ find . | grep "R$" | xargs grep -n "ScaleData" --color=auto


()
泛型函数:

# ./seurat-4.1.0/R/generics.R:565:ScaleData <- function(object, ...) {
#' Scale and center the data.
#'
#' Scales and centers features in the dataset. If variables are provided in vars.to.regress,
#' they are individually regressed against each feature, and the resulting residuals are
#' then scaled and centered.
#'
#' ScaleData now incorporates the functionality of the function formerly known
#' as RegressOut (which regressed out given the effects of provided variables
#' and then scaled the residuals). To make use of the regression functionality,
#' simply pass the variables you want to remove to the vars.to.regress parameter.
#'
#' Setting center to TRUE will center the expression for each feature by subtracting
#' the average expression for that feature. Setting scale to TRUE will scale the
#' expression level for each feature by dividing the centered feature expression
#' levels by their standard deviations if center is TRUE and by their root mean
#' square otherwise.
#'
#' @param object An object
#' @param ... Arguments passed to other methods
#'
#' @rdname ScaleData
#' @export ScaleData
#'
ScaleData <- function(object, ...) {
  UseMethod(generic = 'ScaleData', object = object)
}



()
# ./seurat-4.1.0/R/preprocessing.R:2886:ScaleData.Seurat <- function(



#' @param assay Name of Assay to scale
#'
#' @rdname ScaleData
#' @concept preprocessing
#' @export
#' @method ScaleData Seurat
#'
ScaleData.Seurat <- function(
  object,
  features = NULL,
  assay = NULL,
  vars.to.regress = NULL,
  split.by = NULL,
  model.use = 'linear',
  use.umi = FALSE,
  do.scale = TRUE,
  do.center = TRUE,
  scale.max = 10,
  block.size = 1000,
  min.cells.to.block = 3000,
  verbose = TRUE,
  ...
) {
  # 使用默认 assay，如果没有指定
  assay <- assay %||% DefaultAssay(object = object)
  # 获取 实验对象
  assay.data <- GetAssay(object = object, assay = assay)

  # vars.to.regress 如果有在 meta.data 列名中
  if (any(vars.to.regress %in% colnames(x = object[[]]))) {
  	# 则取出该列数据
    latent.data <- object[[vars.to.regress[vars.to.regress %in% colnames(x = object[[]])]]]
  } else {
  	# 否则，空
    latent.data <- NULL
  }

  # 如果 split.by 是字符串，且长度为1
  if (is.character(x = split.by) && length(x = split.by) == 1) {
  	# 则 split.by 等于 Seurat 对象的[[ 该列数据
    split.by <- object[[split.by]]
  }

  # 传给 Assay 对象的该方法
  assay.data <- ScaleData(
    object = assay.data,
    features = features,
    vars.to.regress = vars.to.regress,
    latent.data = latent.data, #新增这一个变量
    split.by = split.by, #如果输入参数是字符串，长度1，则这里是一列
    model.use = model.use,
    use.umi = use.umi,
    do.scale = do.scale,
    do.center = do.center,
    scale.max = scale.max,
    block.size = block.size,
    min.cells.to.block = min.cells.to.block,
    verbose = verbose,
    ...
  )

  # 覆盖回该 Assay 对象
  object[[assay]] <- assay.data

  # 记录日志
  object <- LogSeuratCommand(object = object)
  return(object)
}







()
# ./seurat-4.1.0/R/preprocessing.R:2828:ScaleData.Assay <- function(

#' @rdname ScaleData
#' @concept preprocessing
#' @export
#' @method ScaleData Assay
#'
ScaleData.Assay <- function(
  object,
  features = NULL,
  vars.to.regress = NULL,
  latent.data = NULL,
  split.by = NULL,
  model.use = 'linear',
  use.umi = FALSE,
  do.scale = TRUE,
  do.center = TRUE,
  scale.max = 10,
  block.size = 1000,
  min.cells.to.block = 3000,
  verbose = TRUE,
  ...
) {
  #(A1) 确定使用哪个 slot
  # 如果不是 linear，值强制改为 T; 否则原始值
  use.umi <- ifelse(test = model.use != 'linear', yes = TRUE, no = use.umi)
  # 如果 use.umi ==T，则 counts，否则 data(默认)
  slot.use <- ifelse(test = use.umi, yes = 'counts', no = 'data')

  # (A1) 使用哪些基因
  # 如果没有指定基因，则使用 高变基因
  features <- features %||% VariableFeatures(object)
  # 如果基因长度为0，使用该slot的全部基因
  if (length(x = features) == 0) {
    features <- rownames(x = GetAssayData(object = object, slot = slot.use))
  }

  # (A3) 设置 Assay 的 slot="scale.data" 的值
  object <- SetAssayData(
    object = object,
    slot = 'scale.data',

    new.data = ScaleData( #又传给 matrix 的同名方法，见下文
      object = GetAssayData(object = object, slot = slot.use), #"matrix" "array" 
      features = features,
      vars.to.regress = vars.to.regress,
      latent.data = latent.data,
      split.by = split.by,
      model.use = model.use,
      use.umi = use.umi,
      do.scale = do.scale,
      do.center = do.center,
      scale.max = scale.max,
      block.size = block.size,
      min.cells.to.block = min.cells.to.block,
      verbose = verbose,
      ...
    )
  )

  # (A4) Parenting() 这个函数干啥的？
  # 修改上一个函数的环境变量，详见 R tips 2,3
  suppressWarnings(expr = Parenting(
    parent.find = "ScaleData.Seurat",
    # 传递的键值对，覆盖父环境中的同名变量
    features = features,
    min.cells.to.block = min.cells.to.block,
    use.umi = use.umi
  ))

  return(object)
}







()
# ./seurat-4.1.0/R/preprocessing.R:2587:ScaleData.default <- function(
这是一个超过 200 行的函数。

#' @importFrom future nbrOfWorkers
#'
#' @param features Vector of features names to scale/center. Default is variable features.
#' @param vars.to.regress Variables to regress out (previously latent.vars in
#' RegressOut). For example, nUMI, or percent.mito.
#' @param latent.data Extra data to regress out, should be cells x latent data
#' @param split.by Name of variable in object metadata or a vector or factor defining
#' grouping of cells. See argument \code{f} in \code{\link[base]{split}} for more details
#' @param model.use Use a linear model or generalized linear model
#' (poisson, negative binomial) for the regression. Options are 'linear'
#' (default), 'poisson', and 'negbinom'
#' @param use.umi Regress on UMI count data. Default is FALSE for linear
#' modeling, but automatically set to TRUE if model.use is 'negbinom' or 'poisson'
#' @param do.scale Whether to scale the data.
#' @param do.center Whether to center the data.
#' @param scale.max Max value to return for scaled data. The default is 10.
#' Setting this can help reduce the effects of features that are only expressed in
#' a very small number of cells. If regressing out latent variables and using a
#' non-linear model, the default is 50.
#' @param block.size Default size for number of features to scale at in a single
#' computation. Increasing block.size may speed up calculations but at an
#' additional memory cost.
#' @param min.cells.to.block If object contains fewer than this number of cells,
#' don't block for scaling calculations.
#' @param verbose Displays a progress bar for scaling procedure
#'
#' @importFrom future.apply future_lapply
#'
#' @rdname ScaleData
#' @concept preprocessing
#' @export
#'
ScaleData.default <- function(
  object,
  features = NULL,
  vars.to.regress = NULL,
  latent.data = NULL,
  split.by = NULL,
  model.use = 'linear',
  use.umi = FALSE,
  do.scale = TRUE,
  do.center = TRUE,
  scale.max = 10,
  block.size = 1000,
  min.cells.to.block = 3000,
  verbose = TRUE,
  ...
) {
  # 检查参数
  CheckDots(...)
  # (A1) 基因名
  # 基因名，用全名
  features <- features %||% rownames(x = object)
  # 用交集
  features <- as.vector(x = intersect(x = features, y = rownames(x = object)))
  
  # (A2) 获取数据子集，按照基因名
  object <- object[features, , drop = FALSE]
  # 获取行名、列名，list形式
  object.names <- dimnames(x = object)
  
  # 取最小值: 该参数，列数
  # 这个变量干啥用的？ //todo
  min.cells.to.block <- min(min.cells.to.block, ncol(x = object))

  # 修改上一层函数的环境内的变量，详见 R tips 2-3
  suppressWarnings(expr = Parenting(
    parent.find = "ScaleData.Assay",
    features = features,
    min.cells.to.block = min.cells.to.block
  ))

  # (A3) split.by 如果为空，则取值T，也就是至少是1组cell
  split.by <- split.by %||% TRUE
  # 把细胞拆分
  split.cells <- split(x = colnames(x = object), f = split.by)

  # 释放内存
  CheckGC()


  # (A4) 如果有要 回归的变量
  if (!is.null(x = vars.to.regress)) {
  	# (B1) 如果 latent.data 为空
    if (is.null(x = latent.data)) {
      # 则建立空数据框，行名 为 数据框的列名 cell id
      latent.data <- data.frame(row.names = colnames(x = object))
    # 如果 非空
    } else {
      # 则按数据框的列名 (colnames) 取 latent.data 的行
      latent.data <- latent.data[colnames(x = object), , drop = FALSE]
      # 重命名 行名
      rownames(x = latent.data) <- colnames(x = object)
    }

    # (B2) 如果要回归的变量有在 数据框 行名(symbol) 中
    if (any(vars.to.regress %in% rownames(x = object))) {
      # 则为其 latent.data 添加新列
      latent.data <- cbind(
        latent.data,
        #有且仅有一个是基因时有bug，比如 vars.to.regress = c("nFeature_RNA", "CD79A")
        t(x = object[vars.to.regress[vars.to.regress %in% rownames(x = object)], ]) 
        #t(x = object[vars.to.regress[vars.to.regress %in% rownames(x = object)], , drop=F]) #改正后
        # PR: https://github.com/satijalab/seurat/pull/5763 
      )
    }

    # (B3) 更新 vars.to.regress，只保留能在 meta.data 和 基因名 中能找到的
    # 如果 latent.data 为空，则该回归步骤跳过
    # Currently, RegressOutMatrix will do nothing if latent.data = NULL

    # 没找到的回归变量
    notfound <- setdiff(x = vars.to.regress, y = colnames(x = latent.data))

    # 如果都没有找到，则报错
    if (length(x = notfound) == length(x = vars.to.regress)) {
      stop(
        "None of the requested variables to regress are present in the object.",
        call. = FALSE
      )
    # 如果有没找到的回归变量，则警告
    } else if (length(x = notfound) > 0) {
      warning(
        "Requested variables to regress not in object: ",
        paste(notfound, collapse = ", "),
        call. = FALSE,
        immediate. = TRUE
      )
      vars.to.regress <- colnames(x = latent.data)
    }


    # (B4) 如果需要提示
    if (verbose) {
      message("Regressing out ", paste(vars.to.regress, collapse = ', '))
    }

    # (B5) 对数据框的行(symbol)，分成若干份，每份 block.size =1000 默认
    chunk.points <- ChunkPoints(dsize = nrow(x = object), csize = block.size)

    # 如果使用超过1个核，future_lapply 实现多进程，先跳过 //todo
    if (nbrOfWorkers() > 1) { # TODO: lapply
      # (C1) 根据因子组合创建数据框
      chunks <- expand.grid( #见 源码解析 11
        names(x = split.cells),
        1:ncol(x = chunk.points),
        stringsAsFactors = FALSE
      )

	#chunks=expand.grid(
	#  c("n1", "n2"),
	#  1:3,
	#  stringsAsFactors = FALSE
	#)
    #
	#> chunks
	#  Var1 Var2
	#1   n1    1
	#2   n2    1
	#3   n1    2
	#4   n2    2
	#5   n1    3
	#6   n2    3


      # (C2) 启动多线程
      object <- future_lapply(
        X = 1:nrow(x = chunks), #对 chunks 进行循环
        FUN = function(i) {
          row <- chunks[i, ] #获取一行。数据框的本质也是 列list
          group <- row[[1]] #第一列
          index <- as.numeric(x = row[[2]]) #第二列
          return(RegressOutMatrix( # 这个函数下文有
            data.expr = object[chunk.points[1, index]:chunk.points[2, index], 
            					split.cells[[group]], drop = FALSE], #按行列式
            latent.data = latent.data[split.cells[[group]], , drop = FALSE],
            features.regress = features, #要回归的基因
            model.use = model.use,
            use.umi = use.umi,
            verbose = FALSE
          ))
        }
      )

      # 如果细胞分类 大于1 (一般不走这里)
      if (length(x = split.cells) > 1) {
      	# 把 1:n 按照每份 by 个分组，每个编号前后都不在一组，见接下来的例子
        merge.indices <- lapply(
          X = 1:length(x = split.cells),
          FUN = seq.int,
          to = length(x = object),
          by = length(x = split.cells)
        )
		# lapply(
		#   X = 1:3,
		#   FUN = seq.int,
		#   to = 20,
		#   by = 3
		# )
		# 输出:
		# [[1]]
		# [1]  1  4  7 10 13 16 19
		# 
		# [[2]]
		# [1]  2  5  8 11 14 17 20
		# 
		# [[3]]
		# [1]  3  6  9 12 15 18
		# #
		# 按列行合并，每行一个基因
        object <- lapply(
          X = merge.indices,
          FUN = function(x) {
            return(do.call(what = 'rbind', args = object[x]))
          }
        )
        # 按列合并
        object <- do.call(what = 'cbind', args = object)
      # 如果细胞没有分类 (一般走这里)
      } else {
      	# 对list按行合并为df
        object <- do.call(what = 'rbind', args = object)
      }


    # 如果只使用1个核（默认流程）
    } else {
      object <- lapply(
        X = names(x = split.cells), #对每个细胞分组进行循环，不设置时默认只有一个分组TRUE
        FUN = function(x) {
          # 如果 verbose ==T 且 细胞分组超过1个，则给出msg
          if (verbose && length(x = split.cells) > 1) {
            message("Regressing out variables from split ", x)
          }

          # 调用主函数
          return(RegressOutMatrix(
            data.expr = object[, split.cells[[x]], drop = FALSE], #df取子集: 一个分组的细胞
            latent.data = latent.data[split.cells[[x]], , drop = FALSE], #潜在回归数据 的子集
            features.regress = features, #要回归的基因
            model.use = model.use,
            use.umi = use.umi,
            verbose = verbose
          ))
        }
      )
      # 按列合并
      object <- do.call(what = 'cbind', args = object)
    }

    # 重新设置行名、列名：使用前面获取的数据
    dimnames(x = object) <- object.names
    # 释放内存
    CheckGC()
  }
  # end of 处理回归变量



  # (A5) 如果 verbose 且 ( scale or center)
  if (verbose && (do.scale || do.center)) {
  	# 创建提示消息 1
    msg <- paste(
      na.omit(object = c(
        ifelse(test = do.center, yes = 'centering', no = NA_character_),
        ifelse(test = do.scale, yes = 'scaling', no = NA_character_)
      )),
      collapse = ' and '
    )
    # 创建提示消息 2
    msg <- paste0(
      toupper(x = substr(x = msg, start = 1, stop = 1)), #首字母大写
      substr(x = msg, start = 2, stop = nchar(x = msg)), #取其余字母
      ' data matrix'
    )
    message(msg)
  }


  # (A6) 定义 scale.function 方法
  # 如果是 dg*Matrix 对象，则使用系数矩阵的 RowScale 方法
  if (inherits(x = object, what = c('dgCMatrix', 'dgTMatrix'))) {
    scale.function <- FastSparseRowScale #默认
    # > class( pbmc_small@assays$RNA@data )
	# [1] "dgCMatrix"
  } else {
    object <- as.matrix(x = object)
    scale.function <- FastRowScale
  }



  # (A7) 如果大于1个核
  if (nbrOfWorkers() > 1) {
  	# 分组，详见上文 回归 多线程部分
    blocks <- ChunkPoints(dsize = length(x = features), csize = block.size)
    chunks <- expand.grid(
      names(x = split.cells),
      1:ncol(x = blocks),
      stringsAsFactors = FALSE
    )

    scaled.data <- future_lapply(
      X = 1:nrow(x = chunks),
      FUN = function(index) {
        row <- chunks[index, ]
        group <- row[[1]]
        block <- as.vector(x = blocks[, as.numeric(x = row[[2]])])
        # 组合参数列表
        arg.list <- list(
          mat = object[features[block[1]:block[2]], split.cells[[group]], drop = FALSE],
          scale = do.scale,
          center = do.center,
          scale_max = scale.max,
          display_progress = FALSE
        )
        arg.list <- arg.list[intersect(x = names(x = arg.list), y = names(x = formals(fun = scale.function)))]
        # 调用函数，传参
        data.scale <- do.call(what = scale.function, args = arg.list)

        dimnames(x = data.scale) <- dimnames(x = object[features[block[1]:block[2]], split.cells[[group]]])
        suppressWarnings(expr = data.scale[is.na(x = data.scale)] <- 0) #把NA都替换为0
        CheckGC()
        return(data.scale)
      }
    )
    #汇总结果
    if (length(x = split.cells) > 1) {
      merge.indices <- lapply(
        X = 1:length(x = split.cells),
        FUN = seq.int,
        to = length(x = scaled.data),
        by = length(x = split.cells)
      )
      scaled.data <- lapply(
        X = merge.indices,
        FUN = function(x) {
          return(suppressWarnings(expr = do.call(what = 'rbind', args = scaled.data[x])))
        }
      )
      scaled.data <- suppressWarnings(expr = do.call(what = 'cbind', args = scaled.data))
    # 默认流程走这个: 只有一个cell分组。按行合并
    } else {
      suppressWarnings(expr = scaled.data <- do.call(what = 'rbind', args = scaled.data))
    }
  # 如果只有一个核(默认)
  } else {
  	# 先造一个填充NA的矩阵，行列数、行名、列名都和输入矩阵一致
    scaled.data <- matrix(
      data = NA_real_,
      nrow = nrow(x = object),
      ncol = ncol(x = object),
      dimnames = object.names
    )
    # 对基因分组，最大分组数
    max.block <- ceiling(x = length(x = features) / block.size)
    # 细胞分组，默认只有1个
    for (x in names(x = split.cells)) {
      if (verbose) {
        if (length(x = split.cells) > 1 && (do.scale || do.center)) {
          message(gsub(pattern = 'matrix', replacement = 'from split ', x = msg), x)
        }
        pb <- txtProgressBar(min = 0, max = max.block, style = 3, file = stderr()) #进度条 step1
      }

      # 对基因的每一个分块(默认1k个基因)
      for (i in 1:max.block) {
      	# 一个分块到另一个分块，整体+1，表示从下标1开始
        my.inds <- ((block.size * (i - 1)):(block.size * i - 1)) + 1
        # 最大不能超过基因个数
        my.inds <- my.inds[my.inds <= length(x = features)]

        # 创建参数列表
        arg.list <- list(
          mat = object[features[my.inds], split.cells[[x]], drop = FALSE],
          scale = do.scale,
          center = do.center,
          scale_max = scale.max,
          display_progress = FALSE
        )

        # 取其与 函数参数 的交集
        arg.list <- arg.list[intersect(
        	x = names(x = arg.list), 
        	y = names(x = formals(fun = scale.function))
        	)]
        # 执行函数
        data.scale <- do.call(what = scale.function, args = arg.list)

        # 行列名字
        dimnames(x = data.scale) <- dimnames(x = object[features[my.inds], split.cells[[x]]])

        # 填充到NA矩阵中
        scaled.data[features[my.inds], split.cells[[x]]] <- data.scale

        # 删除掉中间变量，并释放内存
        rm(data.scale)
        CheckGC()
        
        if (verbose) {
          setTxtProgressBar(pb = pb, value = i) # 进度条 step2
        }
      }

      if (verbose) {
        close(con = pb) #进度条 step3
      }
    }
  }

  # 整体又覆盖一次行列名字 //是不是重复? //todo
  dimnames(x = scaled.data) <- object.names
  # 把矩阵中的NA改为0
  scaled.data[is.na(x = scaled.data)] <- 0
  CheckGC()
  # 返回该矩阵
  return(scaled.data)
}















() Parenting()

# ./seurat-4.1.0/R/utilities.R:2232:Parenting <- function(parent.find = 'Seurat', ...) {
# 注释说：Parenting parameters from one environment to the next
# 修改父环境中的变量，方便日志记录。

# Parenting parameters from one environment to the next
#
# This function allows one to modify a parameter in a parent environment
# The primary use of this is to ensure logging functions store correct parameters
# if they've been modified by a child function or method
#
# @param parent.find Regex pattern of name of parent function call to modify.
# For example, this can be the class name for a method that was dispatched previously
# @param ... Parameter names and values to parent; both name and value must be supplied
# in the standard \code{name = value} format; any number of name/value pairs can be specified
#
# @return No return, modifies parent environment directly
#
# @examples
# Parenting(parent.find = 'Seurat', features = features[features > 7])
#
Parenting <- function(parent.find = 'Seurat', ...) {
  # 获取函数调用栈，最先调用的在最右侧
  calls <- as.character(x = sys.calls())
  # 使用"("分隔，获取第一个元素，也就是调用的函数名历史，第一个是最新的
  calls <- lapply(
    X = strsplit(x = calls, split = '(', fixed = TRUE),
    FUN = '[',
    1
  )

  # 匹配给定的关键词的下标
  parent.index <- grep(pattern = parent.find, x = calls)

  # 如果匹配下标长度不是1，也可能是0 或者2？
  if (length(x = parent.index) != 1) {
  	# 警告：找不到 父环境
    warning(
      "Cannot find a parent environment called ",
      parent.find,
      immediate. = TRUE,
      call. = FALSE
    )
  # 如果下标长度是1
  } else {

  	# 获取传递给父对象的参数列表
    to.parent <- list(...)
    # 如果长度为0，警告：没
    if (length(x = to.parent) == 0) {
      warning("Nothing to parent", immediate. = TRUE, call. = FALSE)
    # 长度不是0时
    # 如果没有name，报错
    } else if (is.null(x = names(x = to.parent))) {
      stop("All input must be in a key = value pair")
    # 如果过滤掉长度为0的names后，和原来个数有差别，报错：所有参数必须有名字
    } else if (length(x = Filter(f = nchar, x = names(x = to.parent))) != length(x = to.parent)) {
      stop("All inputs must be named")
    # 开始干活
    } else {
      # 获取上一个函数栈的环境
      parent.environ <- sys.frame(which = parent.index)
      # 把值逐个写入
      for (i in 1:length(x = to.parent)) {
        parent.environ[[names(x = to.parent)[i]]] <- to.parent[[i]]
      }
    }
  }
}










() RegressOutMatrix()
./seurat-4.1.0/R/preprocessing.R:3366:RegressOutMatrix <- function(


# 回归掉技术效应 和 细胞周期。
# Regress out techincal effects and cell cycle from a matrix
#
# Remove unwanted effects from a matrix
#
# @parm data.expr An expression matrix to regress the effects of latent.data out
# of should be the complete expression matrix in genes x cells
# @param latent.data A matrix or data.frame of latent variables, should be cells
# x latent variables, the colnames should be the variables to regress
# @param features.regress An integer vector representing the indices of the
# genes to run regression on
# @param model.use Model to use, one of 'linear', 'poisson', or 'negbinom'; pass
# NULL to simply return data.expr
# @param use.umi Regress on UMI count data
# @param verbose Display a progress bar
#
#' @importFrom stats as.formula lm
#' @importFrom utils txtProgressBar setTxtProgressBar
#
RegressOutMatrix <- function(
  data.expr,
  latent.data = NULL,
  features.regress = NULL,
  model.use = NULL, #默认使用线性模型 'linear'
  use.umi = FALSE,
  verbose = TRUE
) {
  # (A1) 如果 latent.data 或 模型 为空，则直接返回原始数据框
  # Do we bypass regression and simply return data.expr?
  bypass <- vapply(
    X = list(latent.data, model.use),
    FUN = is.null,
    FUN.VALUE = logical(length = 1L)
  )
  if (any(bypass)) {
    return(data.expr)
  }

  # (A2) 如果不是这3个模型，则报错
  # Check model.use
  possible.models <- c("linear", "poisson", "negbinom")
  if (!model.use %in% possible.models) {
    stop(paste(
      model.use,
      "is not a valid model. Please use one the following:",
      paste0(possible.models, collapse = ", ")
    ))
  }

  # (A3) 如果要回归的基因为空，则使用全部基因 下标
  # Check features.regress
  if (is.null(x = features.regress)) {
    features.regress <- 1:nrow(x = data.expr)
  }
  # 如果是字符串 基因
  if (is.character(x = features.regress)) {
  	# 求交集：要回归的基因、全部基因
    features.regress <- intersect(x = features.regress, y = rownames(x = data.expr))
    # 如果交集为空，报错
    if (length(x = features.regress) == 0) {
      stop("Cannot use features that are beyond the scope of data.expr")
    }
  # 否则，应该是数字。如果数字最大值超过基因个数，报错
  } else if (max(features.regress) > nrow(x = data.expr)) {
    stop("Cannot use features that are beyond the scope of data.expr")
  }

  # (A4) 检查维度：细胞数不相等则报错
  # Check dataset dimensions
  if (nrow(x = latent.data) != ncol(x = data.expr)) {
    stop("Uneven number of cells between latent data and expression data")
  }

  # (A5) 如果 model.use 不是 linear，则 use.umi强制等于T，否则使用原值(默认)
  # 也就是 linear 使用 data，其他使用 counts
  use.umi <- ifelse(test = model.use != 'linear', yes = TRUE, no = use.umi)

  # (A6) 创建用于回归的 方程式
  # Create formula for regression
  vars.to.regress <- colnames(x = latent.data) #要回归掉的列名
  # 和传入参数 features.regress 啥关系？
  # vars.to.regress 是要去掉的批次效应等噪音
  # features.regress 是要 scale 的基因: 从前面看，默认是 HVG；从本函数看，默认是全部基因。

  # 创建 方程式，先连接成字符串，然后强转
  fmla <- paste('GENE ~', paste(vars.to.regress, collapse = '+'))
  fmla <- as.formula(object = fmla)

  # (A7) 如果模型是线性的
  if (model.use == "linear") {
    # In this code, we'll repeatedly regress different Y against the same X
    # (latent.data) in order to calculate residuals.  Rather that repeatedly
    # call lm to do this, we'll avoid recalculating the QR decomposition for the
    # latent.data matrix each time by reusing it after calculating it once

    # 这种情况，我们将重复回归不同的Y 对于相同的X(latent.data)，来计算残差。
    # 相对于 重复调用 lm 来这样做，我们更愿意避免对 latent.data 矩阵反复QR分解，而是计算一次并复用
    # 要点: 只做一次QR分解，然后以后反复使用
    regression.mat <- cbind(latent.data, data.expr[1,]) #新增一列: 表达矩阵的第一行，就一个基因？
    colnames(regression.mat) <- c(colnames(x = latent.data), "GENE")
    qr <- lm(fmla, data = regression.mat, qr = TRUE)$qr
    rm(regression.mat)
  }

  # (A8) 创建空残差矩阵，内容是NA
  # Make results matrix
  data.resid <- matrix(
    nrow = nrow(x = data.expr),
    ncol = ncol(x = data.expr)
  )

  # (A9) 开启进度条
  if (verbose) {
    pb <- txtProgressBar(char = '=', style = 3, file = stderr())
  }

  # (A10) 对每个要 scale 的基因名
  for (i in 1:length(x = features.regress)) {
  	# (B1) 基因名
    x <- features.regress[i]
    # (B2) 拼凑 线性代数 定义的 矩阵: AX=Y，解方程需要矩阵 [AY]，把基因列加到最后
    regression.mat <- cbind(latent.data, data.expr[x, ])
    # (B3) 最后一列的列名 'GENE'
    colnames(x = regression.mat) <- c(vars.to.regress, 'GENE')
    # (B4) 回归后的残差向量
    regression.mat <- switch(
      EXPR = model.use,
      
      # 默认是 线性 模型:
      'linear' = qr.resid(qr = qr, y = data.expr[x,]),

      'poisson' = residuals(object = glm(
        formula = fmla,
        family = 'poisson',
        data = regression.mat),
        type = 'pearson'
      ),

      'negbinom' = NBResiduals(
        fmla = fmla,
        regression.mat = regression.mat,
        gene = x
      )
    )
    # (B5) 把残差放回到原空 残差矩阵中第i行(symbol)
    data.resid[i, ] <- regression.mat

    if (verbose) {
      setTxtProgressBar(pb = pb, value = i / length(x = features.regress)) #进度条
    }
  }

  # (A11) 关闭 进度条
  if (verbose) {
    close(con = pb)
  }

  # (A12) 如果使用的是 UMI
  if (use.umi) { #默认不走这里
    data.resid <- log1p(x = Sweep( #残差矩阵每行 减去该行的 最小值
      x = data.resid,
      MARGIN = 1,
      STATS = apply(X = data.resid, MARGIN = 1, FUN = min), #残差矩阵每行的最小值
      FUN = '-'
    ))
  }

  # (A13) 指定 行列 名
  dimnames(x = data.resid) <- dimnames(x = data.expr)

  # (A14) 返回残差矩阵
  return(data.resid)
}








() FastRowScale() 更快
./seurat-4.1.0/R/utilities.R:758:FastRowScale <- function(

#' Scale and/or center matrix rowwise
#'
#' Performs row scaling and/or centering. Equivalent to using t(scale(t(mat)))
#' in R except in the case of NA values.
#'
#' @param mat A matrix
#' @param center a logical value indicating whether to center the rows
#' @param scale a logical value indicating whether to scale the rows
#' @param scale_max clip all values greater than scale_max to scale_max. Don't
#' clip if Inf.
#' @return Returns the center/scaled matrix
#'
#' @importFrom matrixStats rowMeans2 rowSds rowSums2
#'
#' @export
#' @concept utilities
#'
FastRowScale <- function(
  mat,
  center = TRUE,
  scale = TRUE,
  scale_max = 10
) {
  # inspired by https://www.r-bloggers.com/a-faster-scale-function/
  # 每行的均值
  if (center) {
    rm <- rowMeans2(x = mat, na.rm = TRUE)
  }

  # 求每一行的方差
  if (scale) {
    if (center) {
      rsd <- rowSds(mat, center = rm)
    } else {
      rsd <- sqrt(x = rowSums2(x = mat^2)/(ncol(x = mat) - 1))
    }
  }

  # 执行 center: 就是减去 row均值
  if (center) {
    mat <- mat - rm
  }

  # 执行 scale：除以 row 方差
  if (scale) {
    mat <- mat / rsd
  }

  # 最大值削平
  if (scale_max != Inf) {
    mat[mat > scale_max] <- scale_max
  }
  return(mat)
}








() Sweep() 更兼容
./seurat-4.1.0/R/utilities.R:2382:Sweep <- function(x, MARGIN, STATS, FUN = '-', check.margin = TRUE, ...) {


作者说该函数增强了 base::sweep ，兼容了 R3.x 和 R4.x。

# Sweep out array summaries
#
# Reimplmentation of \code{\link[base]{sweep}} to maintain compatability with
# both R 3.X and 4.X
#
# @inheritParams base::sweep
# @param x an array.
#
# @seealso \code{\link[base]{sweep}}
#
Sweep <- function(x, MARGIN, STATS, FUN = '-', check.margin = TRUE, ...) {
  # 感觉多余，因为区别仅仅是 base::sweep()第一个参数的大小写问题。
  if (any(grepl(pattern = 'X', x = names(x = formals(fun = sweep))))) {
    return(sweep(
      X = x,
      MARGIN = MARGIN,
      STATS = STATS,
      FUN = FUN,
      check.margin = check.margin,
      ...
    ))
  } else {
    return(sweep(
      x = x, #R4.0 是小写
      MARGIN = MARGIN,
      STATS = STATS,
      FUN = FUN,
      check.margin = check.margin,
      ...
    ))
  }
}



() CheckGC() 释放内存
源码解析 14 讲过，不写了。
./seurat-4.1.0/R/utilities.R:1848:CheckGC <- function() {
./seurat-object-4.0.4/R/utils.R:179:CheckGC <- function(option = 'SeuratObject.memsafe') {

# Call gc() to perform garbage collection
#
CheckGC <- function() {
  if (getOption(x = "Seurat.memsafe")) {
    gc(verbose = FALSE)
  }
}


